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METHOD FOR PROCESSING SEISMIC DATA TO IDENTIFY 
ANOMALOUS ABSORPTION ZONES 

Cross-reference to related applications 

Not applicable. 

Statement regarding federally sponsored research or development 

Not applicable. 

Background of the Invention 

Field of the Invention 

[00011 The invention relates generally to the field of seismic data processing. More 
specifically, the invention reiaics lu liicuiods for ideutifying economically useful earth 
formations directly from processed seismic data. 

Background Art 

[0002] Seismic exploration techniques are used to locate subsurface earth formations that 
^e likely to produce economically useful materials such as petroleum. Seismic 
exploration techniques include deploying one or more seismic energy sources near the 
earth's surface and deploying an array of seismic sensors at or near tfie surface in tiie 
vicinity of the one or more seismic sources. Seismic energy propagates downwardly 
from the source, where it may be reflected by subsurface acoustic impedance boundaries. 
The reflected seismic energy is detected by the sensors in the array. The sensors generate 
electrical and/or optical signals corresponding to the detected seismic energy. The 
signals are typically recorded for processing. 

[0003] Seismic processing known in the art includes determining structures of the 

subsurface earth formations. Typically, structures are inferred by analysis of the two-way 
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travel time of the seismic energy from the source to the various reflective boundaries 
beneath the surface and back to the sensors at or near the surface. 

0004] It is also known in the art to determine various petrophysical properties of the 
subsurface earth formations by analysis of the frequency content of the detected seismic 
energy and the phase and amplitude relationships between the seismic energy generated 
by the source and the seismic energy detected by the sensors. Such analysis includes 
determining one or more seismic "attributes" of the earth formations. Attributes may be 
computed prestack or poststack. Prestack means processing prior to summing or 
"stacking" individual sensor recordings ("traces") according to a predetermined 
relationship, such as common mid point (CMP) or common depth point (CDP). 
Poststack refers to processing after individual sensor recordings have been siimmed or 
stacked. Poststack attributes include, for example, reflection intensity, instantaneous 
frequency, reflection heterogeneity, acoustic impedance, velocity, dip, depth and 
azimuth. Prestack attributes include moveout parameters such as amplitude-versus-ofifset 
(AVO), and interval and average velocities. Further, attributes may be categorized as 
either instantaneous attributes, wavelet attributes or geometrical attributes. Instantaneous 
attributes are attributes whose values are obtained for each data point in the seismic data 
or within a small time window of data points (e.g., a few milliseconds), such as 
amplitude, phase, frequency and power. "Data points" within seismic data typically 
refers to numbers each representing a value of seismic signal amplitude at the instant in 
time at which each of the amplitude values is recorded. Wavelet attributes are the 
instantaneous attributes computed at tiie maximum point of the envelope. The physical 
meaning of all the wavelet attributes is essentially the same as their instantaneous 
. attribute counterparts. Geometrical, or interval, attributes are attributes of a seismic trace 
within a seismic interval which are computed from the reflection configuration and 
continuity. The following references describe aspects of seismic attributes and their 
applications. 

[0005] U.S. Pat. No. 5,226,019 issued to Bahorich states that with reference to seismic 
attributes, "combining multiple (i.e. two or more) descriptors through addition, 
subtraction, multiplication and ratio, or other means can also be successfully employed", 
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and suggests the use of "a product of the average instantaneous amplitude and average 
instantaneous frequency." 

[0006] U.S. Pat. No. 5,884,229 issued to Matteucci, discloses a statistical method for 
quantitatively measuring the lateral continuity of the seismic reflection character of any 
specified location in a subsurface target formation. 

[0007] U.S. Pat No. 5,930,730 issued to Marfurt et al., discloses a system for forming a 
seismic attribute display from calculated measures of semblance and conesponding 
estimates of true dip and true dip azimuth of seismic traces within an analysis cell. 

[0008] U.S. Pat, No. 6,012,018 issued to Hombuckle, relates to a system for identifying 
volumetric subterranean regions bounded by a surface in which a specific seismic 
characteristic has a constant value. It is stated in the '018 patent that, "in a geological 
region where physical characteristics (e.g., the presence of oil or gas) are well-correlated 
with seismic attributes, (e.g., seismic amplitude data), the identification of a subvolume 
bounded by a constant-seismic-attribute-value surface may provide a very useful 
predictor of the volumetric extent of the attribute and hence of the characteristic." 

[0009] U.S. Pat. No. 5,001,677 issued to Masters, discloses a system which treats 
measured attributes derived firom seismic data as components of a vector, estimates a 
background vector representing typical background geologic strata, and then calculates a 
new attribute. As stated in the '677 patent, the preferred embodiment combines 
information about P (compressional) and S (shear) impedance contrasts so as to 
discriminate prospective reservoir strata from surrounding non-reservoir or l)ackground 
strata. 

[0010] U.S. Pat. No. 5,724,309 issued to Higgs et al, discloses a system in which two 
new seismic attributes (dip magnitude and dip azimuth) are derived from instantaneous 
phase. The system comprises determining a spatial frequency value by taking the 
directional spatial derivative of the instantaneous phase for each of a plurality ofx, y, t(z) 
data points in the seismic data and posting the spatial frequency values to identify 
changes within the earth's subsurface. 
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[0011] U.S. Pat. No. 5,870,691 issued to Partyka et al, discloses a method for processing 
seismic data to identify thin beds. Although it is generally recognized that specific 
seismic attributes are related to specific subsurface properties, a need continues to exist 
for advancements in the use of seismic attributes to improve the delineation of subsurface 
regions of the earth to assist in the exploration and production of oil, natural gas and 
other minerals. There is continuing interest in methods for analyzing seismic data so as 
to provide direct indication of the presence of petroleum beneath the earth* s surface. 

Summary of the Invention 

[00121 One aspect of the invention is a method for identifying zones anomalously 
absorptive of seismic energy. The method includes joint time-fi^equency decomposing 
seismic traces. The decomposed traces are low frequency bandpass filtered to determine 
a general trend of mean frequency and bandwidth of the seismic traces. The decomposed 
traces are then high frequency bandpass filtered to determine local variations in the mean 
fi-equency and bandwidth of the seismic traces. Anomalous absorption zones are 
determined where there is difference between the general trend and the local variations. 

[0013] In one embodiment of the method, the decomposition includes Gabor-Morlet 
decomposition. In one embodiment, the bandpass filtering includes averaging over a 
selected length time window. In a preferred embodiment, the decomposed traces are 
normalized by matching the envelope of the decomposed traces to the envelope of the 
originally recorded traces. 

[0014] In a specific embodiment, the decomposed traces are spectrally balanced, and a 
set of relative acoustic impedances is determined using inversion processing. A product 
of the anomalous absorption indication and the relative acoustic impedance is used to 
determine hybrid attributes. 

[0015] Another aspect of the invention is computer program stored in a computer 
readable medium. The program includes logic operable to cause a programmable 
computer to perform the following. First, joint time-fi-equency decomposing seismic 
traces is performed. The decomposed traces are then low frequency bandpass filtering to 
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determine a general trend of mean frequency and bandwidth of the seismic traces. The 
decomposed traces are then high frequency bandpass filtered to determine local 
variations in the mean frequency and bandwidth of the seismic traces, whereby 
anomalously al)sorptive zones are determined when the local variations deviate from the 
general trend. 

[0016] Other aspects and advantages of the invention will be apparent from the following 
description and tiie appended claims. 

Brief Description of the Drawings 

[0017] Figure 1 shows a flow chart of example embodiments of a method according to 
the invention. 

[0018] Figure 2 shows an example of a seismic trace, a trace envelope, spectral 
decomposition and trend analysis according to one embodiment of the invention. 

Detailed Description 

[0019] It is known in the art that as seismic energy propagates through earth formations, 
higher frequ^cy components in tiie seismic energy lose tiie ^ame fractional amount of 
energy per cycle as do lower frequency components of the seismic energy. However, the 
higher frequency components lose more energy over the same travel distance as do lower 
frequency seismic energy components. This results in a shift in the spectral content of 
the seismic energy to lower frequencies, as well as narrowing the bandwidth of the 
seismic energy, as the seismic energy propagates through the earth formations. If all 
earth formations were homogenous, the energy loss could be expressed in the form of a 
monotonic exponential relationship. Any variation of the acoustic properties of the 
various formations within the earth will perturb this general relationship. However, in 
general, the overall energy loss trend will approximate a smooth exponential relationship 
due to integration of the absorption effects over the travel distance. Any anomalously 
high (or low) absorption zone will perturb this smooth trend and will appear as an 
increase (or decrease) in the rate of energy loss with respect to the overall energy loss 
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trend. Thus, the frequency content and bandwidth of the seismic energy with respect to 
depth in the earth, or length of propagation path, can be used in various embodiments of a 
method according to the invention to detect anomalous seismic energy absorption zones. 
In embodiments of a method according to the invention, anomalous absorption zones can 
be detected as a deviation from an overall energy absorption trend. 

[0020] Seismic data may be processed in various embodiments of the invention by 
performing the procedures explained below on a number of individual seismic sensor 
data recordings or "traces.'' A trace typically represents a record, with respect to time, of 
seismic signals detected by a single sensor in an array of sensors for a single actuation of 
a seismic energy source. Referring to Figure 1, at 10, user input parameters are selected, 
such as low and high frequency filter passbands for bandpass filters, a number of 
frequency bands, and upper and lower frequencies for the frequency bands for Gabor- 
Morlet decomposition. The low and high frequency bandpass filters will be fiirther 
explained below. At 12, preprocessing, such as static correction, common depth point 
(CDP) stacking, or time migration may be performed on seismic traces to be processed 
according to the invention. An example of a seismic trace is shown at 30 in Figure 2. 

[0021] Returning to Figure 1, at 14, a joint time-frequency decomposition is performed 
on selected seismic traces after preprocessing. Let S(t) represent the amplitude with 
respect to time of a seismic trace, let and E(t) represent its envelope. An example of a 
trace envelope is shown at 32 in Figure 2. The envelope E(t) represents the energy 
amplitude in the trace with respect to time. In the present embodiment the seismic trace 
S(t) is decomposed into Gabor sub-bands using the Gabor-Morlet decomposition. See, 
for example, T. M. Taner, Joint Time/Frequency Analysis, Q Quality Factor and 
Dispersion Computation Using Gabor-Morlet Wavelets or the Gabor-Morlet Transform, 
Rock Solid Images, Houston, TX (2001). The Gabor-Morlet decomposition can be 
determined by the expression: 

[0022] G(/, /) = Y.S(i - r).g(/, r) (1) 

T 

[0023] where 
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[00241 g(/,r) = e-"'^'.e-^^''' (2) 

[0025] Because the output of the decomposition is a complex trace in each sub-band, 
then: 

[0026] a(/,/) = |G(/,0| (3) 

[0027] where a(f, t) represents the joint time-amplitude spectrum of the seismic trace, t 
represents time and / represents frequency. An example plot showing amplitude of each 
frequency component (in each Gabor-Morlet sub-band) with respect to time is shown in 
Figure 2 at 34. The plot is in the familiar "wiggle trace" format, wherein amplitude of 
the component is indicated by the size of the right-hand deflection of each sub band 
amplitude curve and area fill with respect to a selected reference level. 

[0028] At 16, the bandwidth and the RMS frequency of the seismic traces are determined 
for each time sample, and filter passbands for a low frequency (long time window) and a 
hi^ frequency (short time window) filter are determmed. Determining the bandwidth 
and RMS frequency of the decomposed traces in the present embodiment includes 
determining the first and second moments of the joint time-frequency spectrum of each 
trace (determined by the decomposition as explained above). The first moment 
represents the mean value of the instantaneous frequency spectrum and can be calculated 
by the expression: 

[0029] , ,,.5;^ (4) 

[0030] The second moment represents the instantaneous RMS frequency and can be 

calculated by the expression: 

ir.«(/.o 

[0031] U^ = -^ 

[0032] The difference of the squares of the second and first moments (see equation (6) 
below) is equal to the variance of the frequency spectrum of each trace at each time 
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sample. For a Gaussian-distributed spectrum, the variance is equal to the square of half 
the bandwidth of the seismic energy. In general, therefore, the variance is proportional to 
the bandwidth. The variance can be determined by the expression: 

[0033J CTHt) = FLsit)-FLit) (6) 

[0034] The bandwidth and RMS frequency calculations are instantaneous values for each 
time sample point and may contain local noise. In one embodiment, a general seismic 
energy loss trend is estimated using a weighted low pass filter of the variance (instead of 
actual bandwidth) and of the RMS frequency. Since both the bandwidth and RMS 
frequency generally decrease with respect to time, their product can be used rather than 
the individual bandwidth and RMS frequency values. Using the product of two values is 
equivalent to a logical "AND" condition, thus making the low pass filter more effective 
in reducing noise in the trend calculation. The low pass filter is preferably weighted 
using the trace envelope Eft). Weighting the filter using the trace envelope will provide 
more wei^t where there is high energy in each trace and lower weight where there is low 
energy in each trace. Envelope weighting of the filter is based on the assumption that any 
noise is in the background, generally in relatively low amplitude (and thus energy) zones. 
Mid therefore the effects of noise on the low pass filtered frequency and RMS frequenqr 
will be reduced if lower weight is given to low energy zones. 

[0035] Let C(t) represent the product of RMS frequency and the bandwidth: 

[00361 C(0 = 2.F^(0.fT(0 (7) 

[0037] The envelope weighted filtered output, CW(t), which represents the general 
energy loss trend, wherein the filter operator is represented by W(t), is given by the 
expression: 

Y^Eit-T)C{t-T)W{T) 

[0038] cW(t) = -!-=:, 

X 

[0039] In one embodiment, a selected number of seismic traces are averaged prior to 
lows pass filtering to calculate the general trend In one embodiment, in vy^ch seismic 



8 



PATENT APPLICATION 
ATTORNEY DOCKET NO. RSI-03-02 



traces are acquired in three dimensions, traces along five in-line and five cross-line 
receiver steps (total of twenty five individual traces) are averaged prior to low pass 
filtering. An example of a general trend is shown at curve 36 in Figure 2. One 
embodiment of bandpass filtering includes averaging the values of the bandwidth/RMS 
fi-equency product in a time window having a selected length. A typical filter length for 
the low fi-equency bandpass filter is about 800 milliseconds, although the actual value 
used will depend on the length (recorded time) of the traces and the frequency content of 
the traces. The output of the low frequency bandpass filter represents the general trend of 
bandwidth and RMS frequency of the seismic energy with respect to time. 

[0040] Local variations of the RMS frequency and bandwidth are also computed by 
envelope weighted bandpass filtering, but using a higher frequency bandpass filter. Just 
as for the low frequency filter, the high frequency filter may be a time-averaging window, 
but of shorter time duration. A typical filter window length for the higher frequency 
bandpass filter is about equal to the seismic wavelet time. An example of a local trend 
curve is shown at 38 in Figure 2. 

[0041] Returning once again to Figure 1, having bandpass filtered the decomposed 
traces, using both high frequency and low frequency filters, at 18, anomalous absorption 
zones may be determined as explained below. Let CL(t) and CS(t) represent the low- 
frequency and high-frequency bandpass filtered outputs, respectively. Then the 
difference between the general trend (the output of the low frequency bandpass filter) and 
the local variation (the output of the high frequency bandpass filter) is given by the 
expression: 

[0042] AZ{t) = CS{t) - CL(t) (9) 

[0043] Zones having negative values zones of AZ(t) represent areas of lower frequency 
and narrower bandwidth than would be predicted by the general trend. Negative values 
of AZ(t) can be interpreted as an indication of anomalously high absorption zones. Such 
zones have been known to correspond to petroleum-bearing earth formations. An 
example of such as zone is shown at 40 in Figure 2. Conversely, positive values of AZ(t) 
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may be interpreted as anomalously low absorption zones. The values of AZ(t) may be 
output at 20 in Figure 1 as an anomalous absorption zone indicator. 

[0044] In a particular embodiment, hybrid seismic attributes can be determined Such 
hybrid attributes may include anomalously high absorption/low acoustic impedance 
zones; anomalously high absorption/high acoustic impedance zones; anomalously low 
absorption/low acoustic impedance zones; and anomalously low absorption/low acoustic 
impedance zones. To determine the foregoing hybrid seismic attributes, it is necessary to 
estimate the relative acoustic im^pedance of the various earth formations illmninated by 
the seismic energy during trace acquisition. Returning to Figure 1, one embodiment of 
determining relative acoustic impedance from seismic traces includes, at 15, spectrally 
balancing the previously joint decomposed traces. In the present embodiment, spectral 
balancing includes determining for each time sample the envelope of the decomposed 
traces, and then matching the envelope of the decomposed trace to the envelope of the 
original trace at each time sample. Next, shown generally at 22 in Figure 1, a relative 
acoustic impedance series is determined by "band limited inversion'' of the spectrally 
decomposed traces, can be explained as follows. 

[0045] Amplitude versus offset (AVO) computation techniques known in the art provides 
estimates of pressure wave, shear wave and pseudo-Poissori's reflectivity. Ail of these 
estimates are based on the Aki-Richards approximation of Zoepprits formulation of 
reflection amplitude and polarity variation with respect to incidence angle. The Aki- 
Richards approximation assumes that interfaces have gentle impedance contrast. For 
example, the pressure wave impedance contrast (reflectivity) is defined therein as; 

[0046] R^=:Av^p/v^p (10) 

[0047] where; 

[0048] ^"'^ " '^P^^lPn.l - ^p^Pn 

Vp=(Vp^l/?„^l+V^^/7„)/2. 

[0049] It is desired to compute corresponding velocities, given a set of impedances. 

Since all of the impedances have same type of description, here we will omit the p to give 
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a general expression. By substituting equation (11) back into equation (10) provides the 
expression: 

100501 ^.i=2(v„,,p,,i--v„/?J/(v,,,p„,i+v„p„) . (12) 

[0051] It is possible to solve equation (12) recursively. However, recursive solutions can 
be unstable, therefore they require extreme care in selection and control of scale factors 
of the input data. The exact scales, however, are not known. Equation (12) can be 
modified to include the unknown scale factor c: 

[00521 ^.^.1 == 2iv„,,p„,, - v„p„)/iv„^,p„,, + v„p„) (13) 

[00531 The subsequent layer velocity can be solved from the previous layer velocity as: 

[00541 c.i?„^,(v„,ip„,i + v„pj = 2(v„,i/7„^i - v„pj (14) 

[00551 Combining like components provides the expression: 

[00561 (O'K^i - 2)v^ip„^i = -(c./?„^i + 2)v„p„ (15) 

[00571 Finally, the recursive solution can be expressed as: 

[00581 v„,,p„,, = v„p„.(2 + c.i^,,)/(2-c.i?„,,) (16) 

[00591 It is necessary to determine the constant c that will keep the recursion stable. 
Computed results will be similar to the relative acoustic impedance, therefore the solution 
will represent medium length variation of the physical quantities (the actual acoustic 
impedances). Therefore, the results can be scaled according to their intended use. First, 
consider the scalar c that multiplies each output during the recursion: 

[00601 r„,,=(2^c.R^,)/(2-c,R,J (17) 

[00611 This scalar is computed for each data sample and multiplies the present output 
value to estimate the next output value. The output values must be subject to some 
reasonable degree of variation. Equation (17) shows that c,R must be less than 2. Any 
value close to 2, however, will produce a very large scalar that will result in a very much 
larger result for the subsequent output value. This will lead to quick failure of the 
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recursion. It is also the case that the impedance value (the product of velocity and 
density) is being computed, therefore the calculations must result in positive numbers. 
Therefore the r values must be positive. Also, the scalar r cannot be continuously less 
than 1.0, or more than 1.0. If r is continuously less than 1.0, then v will exponentially 
decay to very small values. If r is continuously larger than 1.0, the results will increase 
exponentially and instability will result in the recursion. Therefore, the r scalars will 
continuously vary around a mean value of 1.0. This suggests that a running mean of r 
can be set to 1.0 and the maximum deviation c.R is kept to less than 2.0 (more probably 
less than 1.0). It is known that in most cases RMS values of reflection coeflFicients are in 
a range of about 0.1 to 0.2. Values of reflection coefficient close to unity typically only 
occur at the surface of the water in marine seismic data, and even then only if water 
surface is substantially still. These observations provide a general set of limits for the 
scaling. Note that the acoustic impedances are being computed from seismic data that 
contain only band limited information. Therefore, very long wave and short wave 
impedances would have to be supplied by other measurements, such as well logs and 
normal moveout (NMO) related velocity measurements. However, even band limited, 
inversion does provide reasonable estimates of relative acoustic impedance of the various 
layers of earth formations. This is used in the present embodiment to enable determining 
hybrid attributes as will be explained below. 

[0062] In the recursion procedure, the largest values of the scalars which produce stable 
results are selected. Then, any long wavelength trend is removed, typically with relative 
acoustic impedance computation. 



[0063] From Fatti's approximation: 

I 



[0064] i?(^)=2 



(1 + tan-^ ^)-4 



s 



sin'^, (19) 



[0065] where d represents the incidence angle, the quantities Mp I Ip (the pressure wave 
impedance reflectivity) and AI^/Is (the shear wave impedance reflectivity) are 
estimated. Inversion of these reflectivity values will provide pressure and shear wave 
impedance estimates. If nearly constant or veiy slowly varying density is assumed, the 
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inversion results can provide pressure and shear wave velocity estimates. Since these are 
estimated from band limited seismic data, velocity estimates will also be band limited. 
Longer wavelength portions may be supplemented by the NMO or migration velocities 
computed from the travel times. See, for example, Fatti, J. L., Vail, P. J., Smith, G. C, 
Strauss, P. J. and Levitt, P. R., 1994, Detection of gas in sandstone reservoirs mingAVO 
analysis: A 3-D seismic case history using the geostack technique. Geophysics, 59, no. 
09, 1362-1376. 

[0066] By definition, pressure wave velocity is given by the expression: 



[0067] = ^(A + 2ju)/p (20) 

[0068] and the shear wave velocity is given by the expression: 



[0069] v^^4juIp (21) 

[0070] where A represents incompressibility, p represents density, // represents 
rigidity, 1^ = v^p represents pressure wave impedance, and = v^p represents shear 
wave impedance. 

[0071] By substituting these back into equations 20 and 21, the products of density with 
incompressibility and rigidity is provided as: 

Ap = ll-2I^^ 

[mi] ^ ' (22) 



[0073] An output of the band limited inversion is a series of relative acoustic impedances 
(relative being defined as positive or negative with respect to a mean value). At 26 in 
Figure 1, the hybrid attributes can be calculated as the product of the relative acoustic 
impedance values determined as explained above and the values of ^Z^i^^ determined as 
explained above with respect to 16 and 18 in Figure 1. The resulting products can be 
used to identify anomalously high absorption/low acoustic impedance zones; 
anomalously high absorption/high acoustic impedance zones; anomalously low 
absorption/low acoustic impedance zones; and anomalously low absorption/low acoustic 
impedance zones as previously explained. 
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[0074] Embodiments of a method according to the invention can provide a way to 
directly infer the presence of hydrocarbon productive zones from seismic data. 
Embodiments of a method according to the invention can also provide hybrid seismic 
attributes useful in the interpretation of seismic data. 

[0075] Another aspect of the invention is related to computer programs stored in 
computer readable media. One embodiment of a program according to this aspect of the 
invention includes logic operable to cause a programmable computer to perform the 

v./iwixiwiii.3 v/J. a iiiwLiiv/11 cio v/yvpiaiiiwi. auwv^ vvitii iwopwi xu x i^uiWiS i ci.ii«j. ^, x. iiv v/v/xxxpux^x 

readable medium may be, without limitation, a floppy disk, CD-ROM, hard drive, or any 
other computer readable storage device known in the art. Any general purpose 
programmable computer, such as a personal computer or work station may be 
programmed to implement methods according to the invention and as described above in 
the example embodiment of Figure 1. 

I W / A/J VTXXXX^ XXXW XXXWXXXXV/XX XX(X>9 XyVS/XX \XWOVXXl/Wl VVXtXX X^OpWVX XV/ CX XXXXXXXV/VX XXUX.XXL/WX ML 

embodiments, those skilled in the art, having benefit of this disclosure, will appreciate 
that other embodiments can be devised which do not depart from the scope of the 
invention as disclosed herein. Accordingly, the scope of the invention should be limited 
only by the attached claims. 
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